###########################################################
#
# Patrick Tan
#
# Code for machine learning using CMA (Classification of Microarray Analysis)
#
# setwd('/home/tan/ubc/stat540-high-dimensional/epi-cancer'); source('code/ml.R')
#
# Refs
# Shipp, Margaret A et al. (2002) Diffuse large B-cell lymphoma outcome prediction by gene-expression 
# profiling and supervised machine learning.; Nat Med, 8: 68-74 
# http://www.chibi.ubc.ca/Gemma/expressionExperiment/showExpressionExperiment.html?id=223
# http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE9327
# http://www.broadinstitute.org/mpr/lymphoma/
# http://lausanne.isb-sib.ch/~darlene/tp-class.html
#
###########################################################

seed      <- 100  # rand seed
k         <- 200  # minimum number of genes / features to select
train.pct <- 0.75  # fraction of training and testing

# install libraries
# source("http://bioconductor.org/getBioC.R")
# getBioC("CMA")
# getBioC("ROCR")
# getBioC("GEOquery")
# install.packages("e1071") # svm

# load libraries
library(ROCR)
library(CMA)
library(GEOquery)
library(gplots)

# well-studied validated gene markers
dlbcl.genes <- c("LDHA", "HMGB2", "LGALS", "CTSD", "CTSB", "RCH1", "KPNA", "CDC47", "TRAP1", "MIF", "CCNB1", "BCL2A1", "CIP2", "HCK", "TFRC")
fl.genes <- c("ITGA4", "SLC", "CLU", "SLC6A9", "DARC", "TCRBV10S1P", "CD3D", "CD3E", "CD40LG", "TNFRSF25")
literature.genes <- c("EZH2", "CIITA", "FOXL2", "BTG", "MEF2B", "CARD", "FOX", "KLH6", "MLL2", "CREBBP", "CBP", "P3000")
validated.genes <- append(dlbcl.genes, fl.genes)
validated.genes <- append(validated.genes, literature.genes)

# 1. load data and classification
# read gemma file format, filtered data
dat.in.file <- 'input/223_shipp-dlbcl_expmat.data.txt'
des.in.file <- 'input/223_shipp-dlbcl_expdesign.data.txt'
dat <- read.delim2(dat.in.file, comment.char="#")
des <- read.delim2(des.in.file, comment.char="#")

# only keep numbers, put gene names out for now
dat.bak        <- dat
rownames(dat)  <- dat$Probe
dat            <- dat[,5:length(colnames(dat))] # cols 1-4 are metadata
rownames(des)  <- des$Bioassay
des            <- des[colnames(dat),]
all(des$Bioassay == colnames(dat))

# write to output file to convert factors into numbers!
# and read them back in as a clean matrix
write.table(dat, file='output/223_shipp-dlbcl_expmat.num.data.txt', sep="\t")
dat <- read.table(file='output/223_shipp-dlbcl_expmat.num.data.txt')

# clean name 
# from: FSCC19onAffymetrixGeneChipHumanFullLengthArrayHuGeneFL___FSCC19onAffymetrixGeneChipHumanFullLengthArrayHuGeneFL 
# to:   FSCC19
colnames(dat) <- sapply(colnames(dat), function(s) { sub('[a-z]+.*', '', s) })
rownames(des) <- sapply(rownames(des), function(s) { sub('[a-z]+.*', '', s) })
all(rownames(des) == colnames(dat))

# scale rows, center at 0, normalize sd
dat <- t(scale(t(dat)))
pdf('output/shipp-scale.pdf', width=10)
par(oma=c(2,2,2,2))
boxplot(dat, ylab='log2 normalized gene expression', main='Boxplot of normalized data', pars=list(cex.axis=0.5, las=2))
dev.off()

dim(dat)  #4372   77
all(!is.na(dat)) # no NA
summary(des) # DiseaseState: 58 Diffuse_large_B_cell_lymphoma, 19 follicular_lymphoma

# batch-effect analysis
# pdf('output/shipp-cor-heatmap.pdf')
# par(oma=c(2,2,2,2))
# cor.mat <- cor(as.matrix(dat))
# cor.mat[which(cor.mat < -0.5)] <- -0.5
# cor.mat[which(cor.mat > 0.5)] <- 0.5
# breaks <- seq(-1, 1, 0.2)
# heatmap.2(cor.mat, trace='none', keysize=1, density.info='none', symm=T, Rowv=T, Colv=T,
#   main=paste('Sample correlation matrix\n red=DLBCL, blue=FL'), cexRow=0.5, cexCol=0.5,
#   RowSideColors=c('red','blue')[des$DiseaseState])
# dev.off()

# 2. separate training and testing samples
# split train and test data, our data is pretty small so this might be a problem ...
# split the data into training (75% ~ 65) and testing (15% ~ 12)
all.dlbc <- rownames(des[des$DiseaseState == "Diffuse_large_B_cell_lymphoma",])
all.fl <- rownames(des[des$DiseaseState != "Diffuse_large_B_cell_lymphoma",])

# split into train and test sets
set.seed(seed); 
train.dlbc <- sample(all.dlbc, round(train.pct * length(all.dlbc)))
set.seed(seed);
train.fl <- sample(all.fl, round(train.pct * length(all.fl)))
train.names.dat <- append(train.dlbc, train.fl)
length(train.names.dat)
train.dat <- dat[,train.names.dat]
train.des <- des[train.names.dat,]
dim(train.dat)  #4372 X 58 is the training, the rest are tests

# 3. GenerateLearningSets (methods=CV, MCCV, bootstrap)
set.seed(seed)
splits <- GenerateLearningsets(y=train.des$DiseaseState, method="CV", fold=10, strat=T) 
# 48 observations per learning set

# 4. GeneSelection (feature selection, variable selection) on those learning sets, using method=t.test, wilcoxin.test, limma
# pick top genes from train data

# CMA expects the samples to be in rows.
varsel_t.test <- GeneSelection(X=t(train.dat), y=train.des$DiseaseState, learningsets=splits, method="t.test")
varsel_wilcox.test <- GeneSelection(X=t(train.dat), y=train.des$DiseaseState, learningsets=splits, method="wilcox.test")
varsel_limma <- GeneSelection(X=t(train.dat), y=train.des$DiseaseState, learningsets=splits, method="limma")

# look at genes we got
toplist(varsel_t.test)
toplist(varsel_wilcox.test)
toplist(varsel_limma)

# aggregate results from all iterations

# # t.test doesn't have good clustering
# seliter<-numeric()
# for(i in 1:10) seliter<-c(seliter, toplist(varsel_t.test, iter=i, k = k, show=FALSE)$index)
# dat.bak[as.numeric(names(sort(table(seliter), dec=T))), c("GeneSymbol")]
# 
# # wilcox.test is much better, has BCL2A1, but bad top genes clustering
# seliter<-numeric()
# for(i in 1:10) seliter<-c(seliter, toplist(varsel_wilcox.test, iter=i, k = k, show=FALSE)$index)
# dat.bak[rownames(sort(table(seliter), dec=T)), "GeneSymbol"] # summarize, genes found during cross validation
# 
# # limma, the best method that shows good clustering for top 200 and bad clustering for bottom 200
# seliter<-numeric()
# for(i in 1:10) seliter<-c(seliter, toplist(varsel_limma, iter=i, k = k, show=FALSE)$index)
# sort(table(seliter), dec=T) # summarize
# dat.bak[as.numeric(names(sort(table(seliter), dec=T))), c("GeneSymbol")]

# output top k genes
#shipp.topgenes <- toplist(varsel_t.test, k=dim(train.dat)[1])
#shipp.topgenes <- toplist(varsel_wilcox.test, k=dim(train.dat)[1])
shipp.topgenes <- toplist(varsel_limma, k=dim(train.dat)[1])
shipp.topgenes$genes <- dat.bak[shipp.topgenes$index, "GeneSymbol"]
shipp.topgenes$genes <- as.character(sapply(shipp.topgenes$genes, function(x) { sub("\\|.*",'',x) })) # "BCL2A1|BCL2A1" -> "BCL2A1"
write.table(file='output/shipp-topgenes.csv', shipp.topgenes, sep="\t")
top.names <- shipp.topgenes$genes[1:k]
shipp.top.names <- top.names

pdf('output/shipp-toplist-heatmap.pdf')
par(oma=c(2,2,2,2))
breaks <- seq(-1,1,0.2)
labRow <- sapply(shipp.topgenes$genes[1:k], function(x) {if (x %in% validated.genes) x else ''} )
heatmap.2(dat[shipp.topgenes$index[1:k],], breaks=breaks, keysize=1, density.info='none', trace='none', 
  labRow=labRow, ylab='Probes', cexRow=0.5, cexCol=0.5, ColSideColors=c('red','blue')[des$DiseaseState], 
  main=paste('Top',k,'genes\n red=DLBCL, blue=FL'))
labRow <- sapply(rev(shipp.topgenes$genes)[1:k], function(x) {if (x %in% validated.genes) x else ''} )
heatmap.2(dat[rev(shipp.topgenes$index)[1:k],], breaks=breaks, keysize=1, density.info='none', trace='none', 
  labRow=labRow, cexRow=0.3, cexCol=0.5, ColSideColors=c('red','blue')[des$DiseaseState], 
  main=paste('Bottom',k,'genes\n red=DLBCL, blue=FL'))
dev.off()

# 5. Classification (classifier=svmCMA, knnCMA, pknnCMA)
# radial kernel svm
# this can be slow, so commented out for now
# for k=200, use cost=50, gamma=0.00125
#tunesvm <- tune(X=train.dat.feat, y=train.des$DiseaseState, learningsets=splits, classifier = svmCMA, kernel = "radial", grids = list(cost = c(10, 50, 100), gamma = 2^{-2:2}))
#tunesvm <- tune(X=train.dat.feat, y=train.des$DiseaseState, learningsets=splits, classifier = svmCMA, kernel = "radial")
# best(tunesvm)
#plot(tunesvm)
cost <- 50
gamma <- 0.00125
tuninglist <- list(grids = list(cost=cost, gamma=gamma, kernel="radial"))

# custom feature selection, instead of using nbgene and genesel
train.dat.feat <- (t(train.dat))[, toplist(varsel_limma, k=k)$index]  
svmclass       <- classification(X=train.dat.feat, y=train.des$DiseaseState, 
                  learningsets=splits, classifier=svmCMA, tuninglist=tuninglist, probability=T)
svmr           <- join(svmclass)

# 6. Visualize results and evaluate results (using ftable(), evaluation(), roc(), obsinfo())
ftable(svmr) # confusion matrix
# number of missclassifications:  1 
# missclassification rate:  0.017 
# sensitivity: 0.929 
# specificity: 1 
#     predicted
# true  0  1
#    0 44  0
#    1  1 13
obsinfo(evaluation(svmclass, scheme="observationwise")) #index 45 misclassified
evaluation(svmclass, measure="auc")

pdf('output/shipp-train-roc.pdf')
#plot(performance(ROCR::prediction(tres@prob[,2], tres@yhat), "tpr", "fpr"))
# Note that misclassifcation performance can (partly widely) differ from the Area under the ROC (AUC). 
# This is due to the fact that misclassifcation rates are always computed for the threshold 'probability = 0.5'.
roc(svmr) #AUC 1
dev.off()

# 7. Apply optimized parameters to full training set 
# one split, code from class
shipp.m          <- matrix(which(rownames(des) %in% train.names.dat), 1)
full.learningset <- new("learningsets", learnmatrix=shipp.m, method="my own", ntrain=length(train.names.dat), iter=1)

# number of missclassifications:  12 (paper only missed 6 samples)
dat.feat         <- (t(dat))[, shipp.topgenes$index[1:k]]  # custom feature selection, instead of using nbgene and genesel
set.seed(seed) # changes, even after setting seed
testclassif      <- classification(X=dat.feat, y=des$DiseaseState, learningsets=full.learningset, 
                    classifier=svmCMA, probability=T, tuninglist=tuninglist)

# Evaluation
tres<-testclassif[[1]]
# tres <- join(testclassif)  #requires probabilities, 
ftable(tres)
# number of missclassifications:  1 
# missclassification rate:  0.053 
# sensitivity: 0.8 
# specificity: 1 
#     predicted
# true  0  1
#    0 14  0
#    1  1  4
obsinfo(evaluation(testclassif, scheme="observationwise"))
evaluation(testclassif, measure="auc")
pdf('output/shipp-roc.pdf')
roc(tres)  
dev.off()

# sanity check, bottom ROC worse than top ROC (using limma, wilcox doesn't have this property)
dat.feat        <- (t(dat))[, rev(shipp.topgenes$index)[1:k]]  # custom feature selection, instead of using nbgene and genesel
set.seed(seed) # changes, even after setting seed
testclassif     <- classification(X=dat.feat, y=des$DiseaseState, learningsets=full.learningset, 
                   classifier=svmCMA, probability=T, tuninglist=tuninglist) #tuneres=tunesvm, k=5 for pknnCMA, nbgene=100

# Evaluation
tres<-testclassif[[1]]
# tres <- join(testclassif)  #requires probabilities, 
ftable(tres)
evaluation(testclassif, measure="auc")
pdf('output/shipp-roc-bottom.pdf')
roc(tres)  
dev.off()

# # null, this takes a while, disable for now
# null <- c()
# sample.size <- 1000
# for (i in 1:sample.size) {
#   s <- sample(1:dim(dat)[1], k)
#   dat.feat <- (t(dat))[, s]  # custom feature selection, instead of using nbgene and genesel
#   testclassif<-classification(X=dat.feat, y=des$DiseaseState, learningsets=full.learningset, 
#          classifier=svmCMA, probability=T, tuninglist=tuninglist) #tuneres=tunesvm, k=5 for pknnCMA, nbgene=100, used cost=10, gamma=0.25 but it seems that any other gene wor
#   tres<-testclassif[[1]]
#   null <- rbind(null, c(i, evaluation(testclassif, measure="misclassification")@score))
# }
# colnames(null) <- c('sample', 'misclassification')
# write.table(null, file='output/shipp-null.csv', sep="\t")
# length(which(null <= 0.053))/sample.size  # 1/19
# length(which(null <= 0.106))/sample.size  # 2/19
# null <- read.table(file='output/shipp-null.csv')
# hist(null[,"misclassification"])
 
#################################################################
# Compare with Ruiz-Vela et al. 2008 DE genes
#################################################################
shipp.genes.all    <- sapply(dat.bak[,"GeneSymbol"], function(x) { sub("\\|.*",'',x) })

gse9327.de         <- read.table(file='output/DE_genes_DLBL.txt', sep="\t")
gse9327.de.bak     <- gse9327.de
gse9327.de$ID      <- sapply(gse9327.de$ID, function(x) { sub(' \\|.*','',x) } )
gse9327.de         <- gse9327.de[which(gse9327.de$ID != "" & gse9327.de$ID %in% shipp.genes.all), ]

par(oma=c(2,2,2,2))
labRow <- sapply(shipp.topgenes$genes[1:k], function(x) {if (x %in% validated.genes) x else ''} )
heatmap.2(dat[shipp.topgenes$index[1:k],], keysize=1, density.info='none', trace='none', 
  labRow=labRow, ylab='Probes', cexRow=0.5, cexCol=0.5, ColSideColors=c('red','blue')[des$DiseaseState], 
  main=paste('Top',k,'genes\n red=DLBCL, blue=FL'))

# number of missclassifications:  12 (paper only missed 6 samples)
de.filter          <- which(shipp.genes.all %in% gse9327.de$ID)
dat.feat           <- (t(dat))[, de.filter]
set.seed(seed) # changes, even after setting seed
testclassif        <- classification(X=dat.feat, y=des$DiseaseState, learningsets=full.learningset, 
                      classifier=svmCMA, probability=T, tuninglist=tuninglist)

# Evaluation
tres<-testclassif[[1]]
# tres <- join(testclassif)  #requires probabilities, 
ftable(tres)
evaluation(testclassif, measure="auc")
pdf('output/shipp-de-roc.pdf')
roc(tres)
dev.off()
obsinfo(evaluation(testclassif, scheme="observationwise"))  #62, 67
des[obsinfo(evaluation(testclassif, scheme="observationwise"))[1]$misclassification$index,]

# overlap between our limma genes with Ruiz-Vela et al. 2008 genes
genes.overlap <- intersect(unique(gse9327.de$ID), unique(shipp.topgenes$genes[1:k]))
length(genes.overlap)
write.table(genes.overlap, file='output/shipp-ruiz-gene-features.txt')

#########################################################################################
# GSE9327 (Ruiz-Vela et a. 2008) FL VS DLBCL
#########################################################################################
gse9327            <- getGEO(file="input/GSE9327_series_matrix.txt.gz", GSEMatrix=T, AnnotGPL=T)

#36-DLBCL GSM237466-GSM237501    
#33-FL GSM237540-GSM237572

dim(exprs(gse9327))

as.character(colnames(exprs(gse9327)))

dlbcl.cols.ori      <- which(as.character(colnames(exprs(gse9327))) %in% sapply(237466:237501, function(x) {paste('GSM',x,sep="")} ))
fl.cols.ori         <- which(as.character(colnames(exprs(gse9327))) %in% sapply(237540:237572, function(x) {paste('GSM',x,sep="")} ))
wt.cols.ori         <- which(as.character(colnames(exprs(gse9327))) %in% sapply(237414:237418, function(x) {paste('GSM',x,sep="")} ))

length(dlbcl.cols.ori)
length(fl.cols.ori)
length(wt.cols.ori)

# extract dlbcl and fl
samples.sel                  <- append(dlbcl.cols.ori, fl.cols.ori)
gse9327.dlbcl.fl.dat         <- exprs(gse9327)[,samples.sel]
gse9327.dlbcl.fl.dat         <- t(scale(t(gse9327.dlbcl.fl.dat)))
gse9327.dlbcl.fl.des         <- pData(gse9327)[samples.sel, c("title","geo_accession")]
gse9327.dlbcl.fl.des$title   <- sapply(gse9327.dlbcl.fl.des$title, function(x) { sub(" [0-9]+","",x) } ) # remove numbers
gse9327.dlbcl.fl.des$title   <- as.factor(gse9327.dlbcl.fl.des$title)
gse9327.dlbcl.fl.GeneSymbol  <- fData(gse9327)[,"Gene.Symbol"]

# na's
na.cols       <- sort(unique(which(is.na(gse9327.dlbcl.fl.dat), arr.ind=T)[,"col"]))
na.row.counts <- apply(gse9327.dlbcl.fl.dat, 1, function(x) { length(which(is.na(x))) } )
na.col.counts <- apply(gse9327.dlbcl.fl.dat, 2, function(x) { length(which(is.na(x))) } )
#heatmap.2(gse9327.dlbcl.fl.dat[names(na.row.counts > 20), ], trace='none')

ignored.genes <- as.numeric(names(sort(na.row.counts, dec=T)))[1:500]

pdf('output/gse9327-missing.pdf')
par(oma=c(2,2,2,2))
heatmap.2(gse9327.dlbcl.fl.dat[ignored.genes,], ylab='Probes', 
  ColSideColors=c('red', 'blue')[gse9327.dlbcl.fl.des$title], trace='none', cexCol=0.5, labR='', 
  Rowv=F, Colv=F, na.col='black', main='Missing values')
legend('top', fill=c('red','blue'),legend=levels(gse9327.dlbcl.fl.des$title))
dev.off();

# exclude data with most missing values (first 500)
# 5632   69
gse9327.dlbcl.fl.dat.bak        <- gse9327.dlbcl.fl.dat
gse9327.dlbcl.fl.dat            <- gse9327.dlbcl.fl.dat[-ignored.genes,]
gse9327.dlbcl.fl.GeneSymbol     <- gse9327.dlbcl.fl.GeneSymbol[-ignored.genes]

# scale
gse9327.dlbcl.fl.dat            <- t(scale(t(gse9327.dlbcl.fl.dat)))
pdf('output/gse9327-scale.pdf', width=10)
par(oma=c(2,2,2,2))
boxplot(gse9327.dlbcl.fl.dat, ylab='log2 normalized gene expression', 
  main='Boxplot of normalized GSE9327 data', pars=list(cex.axis=0.5, las=2))
dev.off()

# cols changed after removing columns
dlbcl.cols    <- which(as.character(colnames(gse9327.dlbcl.fl.dat)) %in% sapply(237466:237501, function(x) {paste('GSM',x,sep="")} ))
fl.cols       <- which(as.character(colnames(gse9327.dlbcl.fl.dat)) %in% sapply(237540:237572, function(x) {paste('GSM',x,sep="")} ))
wt.cols       <- which(as.character(colnames(gse9327.dlbcl.fl.dat)) %in% sapply(237414:237418, function(x) {paste('GSM',x,sep="")} ))

length(dlbcl.cols)
length(fl.cols)
length(wt.cols)

# split train and test data set
all.dlbc        <- dlbcl.cols
all.fl          <- fl.cols
set.seed(seed); 
train.dlbc      <- sample(all.dlbc, round(train.pct * length(all.dlbc)))
set.seed(seed);
train.fl        <- sample(all.fl, round(train.pct * length(all.fl)))
train.names.dat <- append(train.dlbc, train.fl)
length(train.names.dat) #52 training

# 2. do feature selection
set.seed(seed)
splits                          <- GenerateLearningsets(y=gse9327.dlbcl.fl.des$title[train.names.dat], method="CV", fold=10, strat=T) 
gse9327.dlbcl.fl.varsel         <- GeneSelection(X=t(gse9327.dlbcl.fl.dat[,train.names.dat]), 
                                   y=gse9327.dlbcl.fl.des$title[train.names.dat], learningsets=splits, method="limma")
gse9327.dlbcl.fl.topgenes       <- toplist(gse9327.dlbcl.fl.varsel, k=dim(gse9327.dlbcl.fl.dat)[1])
gse9327.dlbcl.fl.topgenes$genes <- as.character(gse9327.dlbcl.fl.GeneSymbol[gse9327.dlbcl.fl.topgenes$index])

write.table(gse9327.dlbcl.fl.topgenes, file='output/gse9327-topgenes.csv', sep="\t")

pdf('output/gse9327-toplist-heatmap.pdf')
k <- 200
par(oma=c(2,2,2,2))
labRow <- sapply(gse9327.dlbcl.fl.topgenes$genes[1:k], function(x) {if (x %in% validated.genes) x else ''} )
heatmap.2(gse9327.dlbcl.fl.dat[gse9327.dlbcl.fl.topgenes$index[1:k],], trace='none', 
    labRow=labRow, cexRow=0.5, cexCol=0.5, keysize=1, density.info='none',
    ColSideColors=c('red','blue')[gse9327.dlbcl.fl.des$title], main=paste('Top',k,'genes\n red=DLBCL, blue=FL'))
labRow <- sapply(rev(gse9327.dlbcl.fl.topgenes$genes)[1:k], function(x) {if (x %in% validated.genes) x else ''} )
heatmap.2(gse9327.dlbcl.fl.dat[rev(gse9327.dlbcl.fl.topgenes$index)[1:k],], trace='none', 
    labRow=labRow, cexRow=0.5, cexCol=0.5, keysize=1, density.info='none',
    ColSideColors=c('red','blue','white')[gse9327.dlbcl.fl.des$title], main=paste('Bottom',k,'genes\n red=DLBCL, blue=FL'))
dev.off()

# 3. try to use features from shipp-2002 on test data set
# one split, code from class
m <- matrix(train.names.dat, 1)
full.learningset<-new("learningsets", learnmatrix=m, method="my own", ntrain=length(train.names.dat), iter=1)

set.seed(seed) # changes, even after setting seed

# gene selection of GSE9327
dat.feat <- (t(gse9327.dlbcl.fl.dat))[, gse9327.dlbcl.fl.topgenes$index[1:k]]  # custom feature selection, instead of using nbgene and genesel
dat.feat[which(is.na(dat.feat))] <- 0
set.seed(seed) # changes, even after setting seed
testclassif<-classification(X=dat.feat, y=gse9327.dlbcl.fl.des$title, learningsets=full.learningset, 
       classifier=svmCMA, probability=T, tuninglist=tuninglist) #tuneres=tunesvm, k=5 for pknnCMA, nbgene=100, used cost=10, gamma=0.25 but it seems that any other gene wor
tres<-testclassif[[1]]
ftable(tres)
evaluation(testclassif, measure="auc")

# # null, p-value testing, commented out because it takes a while ...
# null <- c()
# sample.size <- 1000
# for (i in 1:sample.size) {
#   s <- sample(1:dim(gse9327.dlbcl.fl.dat)[1], k)
#   dat.feat <- (t(gse9327.dlbcl.fl.dat))[, s]  # custom feature selection, instead of using nbgene and genesel
#   dat.feat[which(is.na(dat.feat))] <- 0
#   testclassif<-classification(X=dat.feat, y=gse9327.dlbcl.fl.des$title, learningsets=full.learningset, 
#          classifier=svmCMA, probability=T, tuninglist=tuninglist) #tuneres=tunesvm, k=5 for pknnCMA, nbgene=100, used cost=10, gamma=0.25 but it seems that any other gene wor
#   tres<-testclassif[[1]]
#   null <- rbind(null, c(i, evaluation(testclassif, measure="misclassification")@score))
# }
# colnames(null) <- c('sample', 'misclassification')
# write.table(null, file='output/gse9327-null.csv', sep="\t")
# length(which(null == 0.0))/sample.size
# null <- read.table(file='output/gse9327-null.csv')
